Pumping by flapping in a viscoelastic fluid 
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In a world without inertia, Purcell's scallop theorem states that in a Newtonian fluid a time- 
reversible motion cannot produce any net force or net flow. Here we consider the extent to which 
the nonlinear rheological behavior of viscoelastic fluids can be exploited to break the constraints 
of the scallop theorem in the context of fluid pumping. By building on previous work focusing on 
force generation, we consider a simple, biologically-inspired geometrical example of a flapper in a 
polymeric (Oldroyd-B) fluid, and calculate asymptotically the time-average net fluid flow produced 
by the reciprocal flapping motion. The net flow occurs at fourth order in the flapping amplitude, and 
suggests the possibility of transporting polymeric fluids using reciprocal motion in simple geometries 
even in the absence of inertia. The induced flow field and pumping performance are characterized and 
optimized analytically. Our results may be useful in the design of micro-pumps handling complex 
fluids. 

PACS numbers: 47.57.-s, 47.15. G-, 47.63. Gd 



I. INTRODUCTION 

Imagining oneself attempting to swim in a pool of viscous honey, it is not hard to anticipate that, because of 
the high fluid viscosity, our usual swimming strategy consisting of imparting momentum to the surrounding fluid 
will not be effective. The world microorganisms inhabit is physically analogous to this situation [1]. As a result, 
microorganisms have evolved strategies which exploit the only physical force available to them - namely fluid drag 
- to propel themselves or generate net fluid transport. The success of these propulsion strategics is vital in many 
biological processes, including bacterial infection, spermatozoa locomotion and reproduction, and ciliary transport [2], 
Advances in micro- and nano- manufacturing technology have also allowed scientists to take inspiration from these 
locomotion strategies and design micropumps [3] and microswimmers [4] . 

The fundamental physics of small-scale locomotion in simple (Newtonian) fluids is well understood [5-7] . In contrast, 
and although most biological fluids are non-Newtonian, many basic questions remain unanswered regarding the 
mechanics of motility in complex fluids. Since they usually include biopolymers, most biological fluids of interest 
display rheological properties common to both fluids (they flow and dissipate energy) and solids (they can store 
energy). Examples include the airway mucus, which acts as a renewable and transportable barrier along the airways 
of the lungs to guard against inhaled particulates and toxic substances [8], as well as cervical mucus, which is important 
for the survival and transport of sperm cells [9]. The influence of viscoclasticity of the fluid on cell locomotion has 
been experimentally quantified by a number of studies [10-16], including the change in the waveform, structure, and 
swimming path of spermatozoa in both synthetic polymer solutions and biological mucus [17]. Gastropod mucus is 
another common non-Newtonian biofiuid, which is useful for adhesive locomotion, and its physical and rheological 
properties have been measured [18-20]. Modeling-wise, different constitutive models have been employed to study 
locomotion in complex fluids (see the short review in Ref. [21]). Among these models, the Oldroyd-B constitutive 
equation is the most popular, both because of its simplicity and the fact that it can be derived exactly from kinetic 
theory by modeling the fluid as a dilute solutions of elastic (polymeric) dumbbells [22-25] . Recent quantitative studies 
have suggested that microorganisms swimming by propagating waves along their flagella have a smaller propulsion 
speed in a polymeric fluid than in a Newtonian fluid [21, 26]. Likewise, a smaller net flow is generated by the 
ciliary transport of a polymeric fluid than a Newtonian fluid. Specifically, Lauga [21] considered the problem with 
a prescribed beating pattern along the flagellum, while Fu and Powers [26] prescribed the internal force distribution 
instead; both studies suggest that viscoelasticity tends to decrease the propulsion speed. 
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In a Newtonian fluid, Purcell's scallop theorem states that swimming and pumping in the absence of inertia can 
only be achieved by motions or body deformations which arc not identical under a time-reversal symmetry (so- 
called "non-reciprocal" motion) [1]. This poses of course an interesting challenge in designing artificial swimmers 
and pumps in simple fluids, which has recently been addressed theoretically and experimentally (see the review in 
Ref. [7]). The question we are addressing in this paper is the extent to which the scallop theorem holds in complex 
fluids. Because polymeric fluids display nonlinear rheological properties such as shear-dependance or normal-stress 
differences [23, 24], in general reciprocal motions are effective in polymeric fluids [27]. New propulsion and transport 
methods can therefore be designed on small scales to specifically take advantage of the intrinsic nonlinearities of the 
fluid. The goal of this paper is to study such a system in the context of fluid pumping with a simplified geometrical 
setup where the pumping performance can be characterized analytically. 

For simple flow geometries, it is not obvious a priori whether a simple oscillatory forcing of a nonlinear fluid leads 
to a net (rectified) flow. For example, for all Oldroyd-like fluids, a sinusoidally-forced Couette flow leads to zero time- 
averaged flow [23]. In previous work [28], we considered a biologically- inspired geometric example of a semi-infinite 
flapper performing reciprocal (sinusoidal) motion in a viscoelastic (Oldroyd-B) fluid in the absence of inertia. We 
showed explicitly that the reciprocal motion generates a net force on the flapper occurring at second order in the 
flapping amplitude, and disappearing in the Newtonian limit as dictated by the scallop theorem. However, there was 
no time-average flow accompanying the net force generation at second order [28]. Here, we report on the discovery 
of a net fluid flow produced by the reciprocal flapping motion in an Oldroyd-B fluid. The net flow transport is 
seen to occur at fourth order in the flapping amplitude, and is due to normal-stress differences. The dependence of 
the pumping performance on the actuation and material parameters is characterized analytically, and the optimal 
pumping rate is determined numerically. Through this example, we therefore demonstrate explicitly the breakdown 
of the scallop theorem in complex fluids in the context of fluid pumping, and suggest the possibility of exploiting 
intrinsic viscoelastic properties of the medium for fluid transport on small scales. 

The geometric setup in this paper is motivated by the motion of cilia-like biological appendages. Cilia are short 
flagella beating collaboratively to produce locomotion or fluid transport [5, 29]. For example, cilia cover the outer 
surface of microorganisms such as Paramecium for self-propulsion. They are also present along our respiratory tracts 
to sweep up dirt and mucus and along the oviduct to transport the ova. Our setup is also relevant to the rigid 
projections attached to the flagellum of Ochromonas, known as mastigonemcs, which protrude from the flagcllum into 
the fluid [5] . As waves propagate along the flagcllum, the mastigonemes are flapped back-and-forth passively through 
the fluid, a process which can lead to a change in the direction of propulsion of the organism [30-32] . 

Our study is related to the phenomenon known as steady (or "acoustic") streaming in the inertial realm [33-46], 
which has a history of almost two centuries after being first discovered by Faraday [33]. Under oscillatory boundary 
conditions, as in the presence of an acoustic wave or the periodic actuation of a solid body in a fluid, migration of fluid 
particles occur in an apparently purely oscillating flow, manifesting the presence of nonlinear inertial terms in the 
equation of motion. This phenomenon occurs in both Newtonian and non-Newtonian fluids [37-46]. In particular, it 
was found that the elasticity of a polymeric fluid can lead to a reversal of the net flow direction [37-40] . As expected 
from the scallop theorem, no steady streaming phenomenon can occur in a Newtonian fluid in the absence of inertia. 
However, as will be shown in this paper, the nonlinear rheological properties of viscoelastic fluids alone can lead to 
steady streaming. In other words, we consider here a steady streaming motion arising purely from the viscoelastic 
effects of the fluid, ignoring any influence of inertia. 

Recently, polymeric solutions have been shown to be useful in constructing microfluidic devices such as flux stabi- 
lizers, flip-flops and rectifiers [47, 48]. By exploiting the nonlinear rheological properties of the fluid and geometrical 
asymmetries in the micro-channel, microfluidic memory and control have been demonstrated without the use of ex- 
ternal electronics and interfaces, opening the possibility of more complex integrated microfluidic circuit and other 
medical applications [47]. In the setup we study here, we do not introduce any geometrical asymmetries and exploit 
solely the non-Newtonian rheological properties of the polymeric fluid for microscopic fluid transport. 

The structure of the paper is the following. In §11, the flapping problem is formulated with the geometrical 
setup, governing equations, nondimcnsionalization and the boundary conditions. In §111, we present the asymptotic 
calculations up to the fourth order (in flapping amplitude), where the time-average flow is obtained. We then 
characterize analytically the net flow in terms of the streamline pattern, directionality and vorticity distribution 
(§IV). Next, we study the optimization of the flow with respect to the Deborah number (§V). Our results are finally 
discussed in §VL 
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FIG. 1: Geometrical setup and notations for the flapping calculation. A semi- infinite plane flaps sinusoidally with small 
amplitude e around an average position at right angle with an infinite wall. 



II. FORMULATION 
A. Geometrical setup 

In this paper, we consider a semi-infinite two-dimensional plane flapping sinusoidally with small amplitude in a 
viscoelastic fluid. The average position of the flapper is situated perpendicularly to a flat wall with its hinge point 
fixed in space (see Fig. 1). The flapper is therefore able to perform reciprocal motion with only one degree of freedom 
by flapping back-and-forth. Such a setup is directly relevant to the unsteady motion of cilia-like biological appendages 
(see §1). 

It is convenient to adopt planar polar coordinates system for this geometrical setup. The instantaneous position of 
the flapper is described by the azimuthal angle 6(t) = ir/2 + e®(t), where Q(t) is an order one oscillatory function of 
time and e is a parameter characterizing the amplitude of the flapping motion. The polar vectors e r (#) and e#((9) are 
functions of the azimuthal angle, and the velocity field u is expressed as u = u r e r +ugeg. In this work, we derive the 
velocity field in the the domain (0 < 6 < n/2) in the asymptotic limit of small flapping amplitude, i.e. e C 1; the 
time-averaged flow in the domain (tt/2 < 9 < n) can then be deduced by symmetry. 



B. Governing equations 

We assume the flow to be incompressible and the Reynolds number of the fluid motion to be small, i.e. we neglect 
any inertial effects. Denoting the pressure field as p and the deviatoric stress tensor as r, the continuity equation and 
Cauchy's equation of motion are respectively 

V-u = 0, (1) 
Vp = V-r. (2) 

We also require constitutive equations, which relate stresses and kinematics of the flow, to close the system of 
equations. For polymeric fluids, the deviatoric stress may be decomposed into two components, r = t s + t p , where 
t s is the Newtonian contribution from the solvent and r p is the polymeric stress contribution. For the Newtonian 
contribution, the constitutive equation is simply given by t s = T)^, where r] s is the solvent contribution to the viscosity 
and 7 = Vu + * Vu. For the polymeric contribution, many models have been proposed to relate the polymeric stress 
to kinematics of the flow [22-25]. We consider here the classical Oldroyd-B model, where the polymeric stress, t p , 
satisfies the upper-convected Maxwell equation 

r p + A ^ = Vp j, (3) 
where rj p is the polymer contribution to the viscosity and A is the polymeric relaxation time. The upper-convected 
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derivative for a tensor A is denned as 

v 8 A 

A=^-+u-VA-( t Vu-A + A-Vu), (4) 

at v ' 

and represents the rate of change of A in the frame of translating and deforming with the fluid. From Eq. (3), we 
can obtain the Oldroyd-B constitutive equation for the total stress, r, as given by 

r + Ai t= 7] + A 2 7^ , (5) 

where r\ = rj s + rj p , Ai = A, and A 2 = r) s \/r]. Here, Ai and A 2 are the relaxation and retardation times of the fluid 
respectively. The relaxation time is the typical decay rate of stress when the fluid is at rest, and the retardation time 
measures the decay rate of residual rate of strain when the fluid is stress-free [23, 24]. It can be noted that A 2 < Ai, 
and both are zero for a Newtonian fluid. 



C. Nondimensionalization 



Periodic flapping motion with angular frequency uj is considered in this paper. Therefore, we nondimensionalize 
shear rates by u> and stresses by r/ui. Lengths are nondimcnsionalized by some arbitrary length scale along the flapper. 
Under these nondimensionalizations, the dimcnsionless equations are given by 

V • u = 0, (6a) 

Vp = V • t, (6b) 

v v 

r + Dei t = 7 + De 2 7, (6c) 

where Dei = Aicj and De 2 = A 2 w are defined as the two Deborah numbers and we have adopted the same symbols 
for convenience. 



D. Boundary conditions 

The boundary condition in this problem can be simply stated; on the flat wall (9 = 0), we have the no-slip and the 
no-penetration boundary conditions. In vector notation, we have therefore 

u(r, 9 = 0) = (7) 

along the flat wall. 

Along the flapper, we also have the no-slip condition, u r (r,9 = ir/2 + eO(i)) = 0. The other boundary condition 
imposed on the fluid along the flapper is given by the rotation of the flapper, ug(r, 9 = it/2 + eO(<)) = rfl(t), where 
£l(t) = e0. In vector notation, we have then 

u{r,9 = ir/2 + e6(t)) =rSl(t)e e . (8) 



III. ANALYSIS 



Noting that a two-dimensional setup is considered, the continuity equation, V • u = 0, is satisfied by introducing 
the streamfunction ^(r, 9) such that u r = (8^/89) jr and ug = —d^ /dr. The instantaneous position of the flapper is 
described by the function 9 = 7r/2 + e6(t), and we consider here a simple reciprocal flapping motion with Q(t) = cost. 
Since small amplitude flapping motion (e <c 1) is considered, we will perform the calculations pcrturbatively in the 
flapping amplitude and seek perturbation expansions of the form 

{u,*,x,p,cr} = 6{ui,*i,Ti,pi,CTi} + e 2 {u 2 ,* 2 ,T 2 ,p 2 ,er 2 } + ..., (9) 

where er = —pi + t is the total stress tensor and all the variables in Eq. (9) are defined in the time-averaged domain 
< 9 < tt/2. Since a domain-perturbation expansion is performed, careful attention has to be paid on the distinction 
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between instantaneous and average geometry. Recall that the polar vectors e r (6(t)) and &g(6(t)) are functions of 
the azimuthal angle which oscillates in time. To distinguish the average geometry, we denote (t) = e r (7r/2) and 
(n) = e§(7r/2) as the average directions along and perpendicular to the flapper axis (See Fig. 1). In this paper, (...) 
denotes time-averaging. 

In addition, we employ Fourier notation to facilitate the calculations. In Fourier notation, the actuation becomes 
= Re{e rf } and O = Ke{ie lt }. Because of the quadratic nonlinearities arising from boundary conditions and the 
constitutive modeling, the velocity field can be Fourier decomposed into the anticipated form 

U! = Rc{uie ?;t }, (10a) 

u 2 =Rc{4 0) + u 2 2) e 2it }, (10b) 

u 3 = Refuse" + uf ) e 38 *}, (10c) 

u 4 = Re{uf + ufe 24 * + ufe 4 *'}, (10d) 

with similar decomposition and notation for all other vector and scalar fields. 

We now proceed to analyze Eq. (6) order by order, up to the fourth order, where the time-average fluid flow occurs. 
The boundary conditions, Eqs. (7) and (8), are also expanded order by order about the average flapper position using 
Taylor expansions. 

A. First-order solution 

1. Governing equation 
The first-order Oldroyd-B relation is given by 



which in Fourier space becomes 



Tl+Dei ^ I + De^, ,11) 

ft - 1±2*S+, (12) 
1 + iDei 

We then note that we have V x V • r = by taking the curl of Eq. (6b). Therefore, we take the divergence and then 
curl of Eq. (12) to eliminate the stress and obtain the equation for the first-order streamfunction 

V 4 *i=0. (13) 

2. Boundary conditions 

At 6 = it/2, the boundary condition at this order is given by 

ui = r0(n), (14) 

which becomes 

Hi = ir(n), (15) 
upon Fourier transformation. We also have the no-slip and no-penetration boundary condition at 9 = 0. 

3. Solution 

The solution satisfying the above equation and boundary conditions is given by 

■ 2 

~ %r 

*i = — (cos 20- 1), (16a) 

IT 

u\ r — — — sin 20, (16b) 

IT 

u xe = — (1- cos 20). (16c) 
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B. Second-order solution 

1. Governing equation 
The second-order Oldroyd-B relation is given by 

i+D 4h-( i+Dc 4h 

= De 2 [ui • - (*Vui • 7i + 7i • Vui)] 

- Dei [ui • Vn - ('Vui ■ ti + n • Vui)] . (17) 
Fourier transforming Eq. (17) and using Eq. (12), we obtain the two harmonics as 

(1 + 2iDei) r 2 2) - (1 + 2iDe 2 ) 7 2 2) 

= Ittj^t [ai • v ^ (tvai ■ ^ + ^ ' v " l)] ' (18) 

and 

= ^T+TrSr t fi i ' V ^ - ^ ' ^ + ^ ' vfi ^ ' (19) 

where the starred variables denote complex conjugates in this paper. Finally, taking the divergence and then curl 
of both Eq. (18) and Eq. (19), and using the knowledge of the first-order solution, we obtain the equation for the 
second-order streamfunctions as simply 

V 4 * 2 2) = 0, (20a) 
V 4 * 2 0) = 0. (20b) 

2. Boundary conditions 
The boundary condition at this order is given by 

u 2 = -e^±-ree<t), (21) 

when evaluated at = tt/2. In Fourier notation and with the first-order solution, the boundary conditions for the 
second-order average flow and the second harmonic read 

-,(o) 



u 2 u; = 0, (22a) 



,-,(2) 



IT 



*2 



-<t). (22b) 



In addition, the no-slip and no-penetration boundary condition are imposed at = 0. 

3. Solution 

The solution satisfying the second-order equation and the boundary conditions is given by 

= 0, (23a) 



2 

~ f2l if ( 1 7T 7T \ 



„~,( 2 ) 



ir 



(cos 20 + - sin 20 - l) , (23c) 

u$j = -j (sin 20 - ~ cos 20 - 20 + |) . (23d) 
As anticipated, there is no time-averaged flow at second order, and we proceed with calculations at higher order. 
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Third-order solution 



Governing equation 



The third-order Oldroyd-B relation is given by 



d 



l+D Cl -| T 3 



1 + De 2 - J73 



: De 2 [ui • V7 2 - (' Vui • 72 + 7 2 ' Vui)] 
Dei [u 2 ■ Vti - (*Vu 2 ■ n + ti • Vu 2 )] 
De 2 [u 2 • V7j - (*Vu 2 • 7i + 7i • Vu 2 )] 
Dei [ui ■ Vr 2 - (*Vui ■ t 2 + t 2 ■ Vui)] . 



(24) 



Aiming at obtaining the average flow at 0(e 4 ), we only need to calculate the first harmonic at 0(e 3 ) since the third 
harmonic will only enter the oscillatory part at 0(e 4 ) (see the fourth-order calculations for details). Therefore, upon 
Fourier transform, we obtain the first harmonic component of Eq. (24) as 



(l+ 2 Dei)T^ 1) -(l + iDe 2 )7^ 1) 



1 Dc 2 - Dei 

2 1 + 2iDci 

1 Dc 2 - Dei 

2 1 - iDei 



uj • V 7 ^ 2) - ( £ Vui ■ 7 2 - 
u 2 ' v 7i - Vu 2 ; • 7i + 7i ' Vu 2 



(2) 



r (2) - * 
7 2 ' 



(25) 



where we have used the constitutive relations given by Eqs. (12) and (18). Taking the divergence and then curl of 
Eq. (25), we obtain the equation for the first harmonic of the third-order streamfunction 



+?,(!) 



V 4 * 



3iDci(Dei - Dc 2 ) cos 29 
r 2 (l - zDei)(l + 2iDei)(l + iDc 2 ) ' 



(26) 



2. Boundary condition 



At 9 = 7r/2, the boundary condition at this order is given by 

„du 2 1 2 9 2 ui 1 ■ 2 
U3 = - 2 6 -Bfi ' 2 r0Q (n) ' 



(27) 



evaluating at 9 = tt/2. The boundary condition for the first harmonic component, in Fourier space, is then given by 

(28) 

At 9 = 0, we also have the no-slip and no-penetration boundary condition. 



u 3 1} -— <t>- T <n>. 



3. Solution 



The solution at this order has the form 



8 



(i) 



u {1) 

U 3,r 



^3.1 



IT 2 71" — 4 \ 7t / i\ 7T 

Y a + —i cos 20 — la + - ) ^29 + - ya , 4 



32 



+ a6» sin 26* 

4-7T 2 



— i - —a] sin 29 - - ( a + - ) cos 29 + - (a + - I + a(29 cos 29 + sin 29) 
16 4/ 2 V 4 / 2 V 4 ' 



7T 7T 

— a+ - 



32 



7T 7T 

-i —a 



a6> sin 26* 



(29a) 
(29b) 



(29c) 



where we have defined the constant 



Q 



3iDci(Dc 2 - Dei) 



8(1 - iDei)(l + 2iDei)(l + iBe 2 ) 



(30) 



D. Fourth-order solution 



Governing equation 



Finally, the fourth-order Oldroyd-B relation is given by 



l + D Cl -jr 4 

= De 2 [ui • V-y 3 

- Dei [ui • Vr 3 - 
+ De 2 [u 2 • V7 2 - (* Vu 2 • 7 2 + 7 2 • Vu 2 ) 

- Dei [u 2 • Vr 2 - (*Vu 2 • r 2 + r 2 • Vu 2 ) 
+ Dc 2 [ua • V 7l - (* Vu 3 • 71 + 7i ■ Vu 3 ) 

- Dei [u 3 • Vri - (*Vu 3 • ti + n • Vu 3 ) 



l + De 2 -j 74 

('VUi-73+73-VU!) 

(*Vui • t 3 + t 3 ■ Vui) 



(31) 



Since we wish to characterize the time-average flow, uf\ we calculate the time-average of Eq. (31) and obtain 



*(<>) _ 

~4 14 



De 2 
Dei 
De 2 



u* • Vx« - 



7 2 2) -V, 



--Dei 
+ ^De 2 
-iDei 



,(2)* 



,(2)* 



-(1)* 
U 3 



Vf 



~(2) 



-,-,(1)* 



,(2) , -(2) 
2 > '2 



4 2) ' 



■V7i-( t V^'.7 1 +7 1 -Vu 3 
• Vfi - f * Vuf • ri + fi • VuW' 



(32) 
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As done previously, we take the divergence and then curl of Eq. (32), and invoke the lower-order constitutive relations 
Eqs. (12), (18) and (25), to obtain the equation for streamfunction of the average flow 

V 4^(0) = p M sin 26 + B 4 cos 26 + C 4 sin 46 ^ 

where 



fl = n„ , r,: 2 ^: (34) 



De 2 - Dei 
2(1 +De?)(2Dei - i)'' 
A 4 = 8a - Dei + 8iaDei + 4iDe^ + 16aDe^, (35) 
B 4 = 27r(l + iDe 1 + 2Dei)(a + Q ! *) , (36) 
C 4 = -2 (8a + 8iaDci + 3iDc^ + 16aDei 

+ 4a* +4iDe x a* + 8De?a*). (37) 

2. Boundary conditions 
At 6 — 7r/2, the boundary condition at this order is written as 

which we then Fourier-transform to obtain the boundary condition for u^ . In addition, since wc arc only interested 



as 

r 



in the time-averaged flow, i.e., real part of the solution Re{u4 }, the boundary condition at 6 = tt/2 can be simplified 

Re{u 4 0) } = - (8 8 ^i Re{a}(t), (39) 
where 

-3Dei(Dei - De 2 )(Dei + De 2 + 2Defoe 2 ) 
B{ai 8(l + Dc2)(l + 4Dc2)(l + Dc2) ' 1 ' 

Finally, as usual, we have the no-slip and no-penetration boundary conditions at 6 = 0. 



3. Solution 



Solving the inhomogeneous biharmonic equation with the boundary conditions above, we obtain our main result, 
namely the analytical formula for the time-averaged flow as 
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FIG. 2: Time-averaged vorticity, (W4), as a function of polar angle 9. (a): Fixed Deborah number (De = 100) and r/ a /rj = 0.1 
(blue solid line), 0.01 (red dashed lines) and 0.001 (black dotted line); (b): Fixed relative viscosity (r/s/ri = 0.1) and De = 1 
(blue solid line), 10 (red dashed line), and 100 (black dotted line). 



Re 



K>} 



Rc 



Rc 



r 2 Dci (Dei - De 2 ) (2De 2 De 2 + Dei + De 2 ) 
512 (De| + 1) 2 (4Dc^ + l) (De 2 . + l) 

32tt - 3tt 3 + 87rDeiDc 2 - 3tt 3 DciDc 2 + 247rDe 2 + 40(-2O + 3tt 2 
- 12Dc 2 - 8DciDc 2 + 37r 2 DeiDe 2 ) + sin 40 (-4 - 12De 2 + 8DeiDe 2 ) 
+ cos 20 (-32tt + 3tt 3 + 486" - 247rDe 2 + 480De 2 - 87rDeiDc 2 + 37r 3 DciDe 2 ) 
+ sin 26» (24 - 6tt 2 + 24De 2 - 247r0De 2 - 67r 2 DeiDe 2 + 247T0DeiDe 2 ) 
rDei (De 2 - Dei) (2De 2 De 2 + Dei + Dc 2 ) 



256 (Def + 1) 2 (4De! 
40 - 6tt 2 

f sin20 (-327r + 37r 3 



1 



(De 2 



1 



24Dc 2 + 16DeiDe 2 - 67r 2 DeiDe 2 + cos 40 1 



24Dc 2 - 16DciDe 2 ) 



480 - 127rDei + 480Dc 2 - 207rDciDc 2 + 37r 3 DeiDe 2 ) 



+ cos 20 (-48 + 6tt 2 - 48De 2 + 247r0De 2 + 67r 2 DeiDe 2 - 247r0DeiDe 2 ) 
rDei (De 2 - Dei) (2De 2 De 2 + Dei + Dc 2 ) 



256 (De 2 + 1) 2 (4De 2 + l) (Dc^ + l) 

32tt - 3tt 3 + 8ttDciDc 2 - 3tt 3 DciDc 2 + 247rDe 2 + 40 (-20 + 3tt 2 

- 12Dc 2 - 8DciDe 2 + 37r 2 DciDe 2 ) + sin 40 (-4 - 12De 2 + 8DeiDe 2 ) 

f cos 20 (-32tt + 3tt 3 + 480 - 24ttDc 2 + 480Dc 2 - 87rDeiDc 2 + 37r 3 DciDe 2 ) 

f sin 20 (24 - 6tt 2 + 24Dc^ - 247r0De 2 - 67r 2 DciDc 2 + 247r0DciDe 2 ) 



(41a) 



(41b) 



(41c) 



IV. CHARACTERIZATION OF THE TIME-AVERAGED FLOW 



In the analysis above, we have computed the flow field perturbatively up to order 0(e 4 ) and found that a nonzero 
time-averaged flow occurs at that order, as described by Eq. (41). Hereafter, for convenience, we rewrite the two 
Deborah numbers as Dei = De and De 2 = De£, where £ = ij s /v is the relative viscosity of the solvent vs. total fluid. 
The creation of a net flow by the tethered flapping motion demonstrates explicitly that Purcell's scallop theorem 
breaks down in a viscoelastic fluid. This suggests that reciprocal flapping-like motion can be exploited for pumping 
polymeric fluids in simple geometries even in the absence of inertia - a situation which is impossible in Newtonian 
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FIG. 3: Streamline and vorticity pattern for r/s/rj = 0.1; (a): De = 1; (b): De = 100. The grayscale map displays the value of 
the vorticity, with legend shown on the right of each plot. 

fluids. In the following sections, we explore the properties of this time-average flow and its dependance on both the 
actuation frequency and material properties of the fluid. 



A. Streamline and vorticity pattern 



With the streamfunction explicitly calculated, we can easily compute the flow streamlines, as well as the flow 
vorticity, given by (w^) = — V (^4), or 



Dc 3 (l - C) (1 + C + 2Dc 2 C) 
128 (1 + De 2 ) 2 (1 + 4De 2 ) (l + De 2 C 2 ) 
-32tt - 24Dc 2 tt + 3tt 3 - 8De 2 7rC + 3De 2 7r 3 C 

- 4 (-20 - 12De 2 + 3tt 2 - 8De 2 C + 3DcVC) 9 

f (24Dc 2 tt - 24De 2 7rC) cos 29 + (48 + 48De 2 ) sin 20 

f- (-12 - 36De 2 + 24Dc 2 C) sin 40 . 



(42) 



We see that the vorticity is only a function of the polar angle 9 between the wall and the flapper. The vorticity 
is plotted as a function of the angle 9 in Fig. 2 for different relative viscosities (Fig. 2a) and different Deborah 
numbers (Fig. 2b). The locations where the vorticity changes its sign are apparently invariant and occur around 
9 w 37r/16 (from negative to positive) and 9 ps 37r/8 (from positive to negative). The streamline pattern and vorticity 
distribution are also qualitatively similar for different Deborah numbers and relative viscosities, as illustrated in Fig. 3 
for different Deborah numbers (De = 1 and De = 100) at a fixed relative viscosity of 0.1. It can be noted that, keeping 
the relative viscosity fixed, increasing the Deborah number leads to more inclined streamlines (greater vertical velocity 
components) near the flat wall (9 = 0). 



B. Directionality of the flow 



As shown from the arrows in the streamline pattern in Fig. 3, the flapping motion draws the polymeric fluid towards 
the hinge point at an acute angle, and pumps the fluid away from the hinge point along both the flat wall and the 
average flapper position. To illustrate this directionality further, we plot in Fig. 4 the radial velocity per unit radius 
against the polar angle for different relative viscosities (Fig. 4a) and Deborah numbers (Fig. 4b). Again, the locations 
where the radial velocity changes its sign are apparently invariant under the change of relative viscosity or Deborah 
number, and occur around 9 ps 7r/4 (positive to negative) and 9 ps 77r/16 (negative to positive). 
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FIG. 4: Net flow velocity along the average flapper position, (ui r )/r, as a function of the polar angle 9. (a): fixed Deborah 
number (De = 100) and r\ a jr\ = 0.1 (blue solid line), 0.01 (red dashed line) and 0.001 (black dotted line); (b) fixed relative 
viscosity {r) a /r) = 0.1) and De = 1 (blue solid line), 10 (red dashed line), and 100 (black dotted line). 



V. OPTIMIZATION 



Having identified the basic flow patterns generated by the flapping motion, we now turn to a possible optimization 
of the pumping performance. Specifically, we address the question: what is the optimal Deborah number at which 
the largest flow can be generated? Since different optimality criteria can be defined, we consider here three different 
"optimality measures" for the net flow, and show they all generate essentially the same conclusion. 



1. Flow along the boundary 

Since the flapping motion pumps the fluid away from the hinge point along the average flapper position, one natural 
measure of the pumping performance is the magnitude of flow along the average flapper position (6 = 7r/2). Note 
that the velocity field is directly proportional to the radius, and recall that the velocity is only radial along the 
average flapper position as required by symmetry. Consequently, the dependence of the intrinsic flow strength upon 
the Deborah number can be characterized by the ratio between the radial velocity along the average flapper position 
and the radial distance, 

(u ir )(r,0 = n/2) 

U b (De, C = 

r 

3Dc 3 (tt 2 -8) (W)(l + C + 2De 2 C) 
64 (1 + 5De 2 + 4De 4 ) (l + De 2 C 2 ) ' 

which is plotted for different relative viscosities in Fig. 5a. From Eq. (43), we see that for small values of De, Ub ~ De 3 , 
whereas for large values of De, Ub ~ 1/De, and therefore an optimal Deborah number is expected to exist. This is 
confirmed in Fig. 5a, where we see that for each value of the relative viscosity, there is an optimal value of the Deborah 
number where the flow along the boundary is maximal. For small relative viscosities, we note the presence of two local 
peaks (in contrast, only one exists for r) s /i] = 1CP 1 ). Physically, by decreasing the relative viscosity, we are varying 
the retardation time of the fluid, while keeping the relaxation time fixed. The position of the second peak changes 
correspondingly and commensurately when the relative viscosity is varied by orders of magnitude, while the position 
the first peak is unchanged. When the relative viscosity is set to zero (zero retardation time, which is a singular 
limit), we see in Fig. 5a a single peak at essentially the same Deborah number as before. From these observations, 
we deduce that the two local optimal Deborah numbers arise from two different properties of the fluid, respectively 
relaxation and retardation. For small relative viscosities, the small local optimal Deborah number can be attributed 
to relaxation while the larger local optimal Deborah number can be attributed to retardation, and it disappears in 
the singular (and unphysical) limit of zero retardation. 
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FIG. 5: Dependence of pumping performance with the Deborah number, for two different pumping measures, (a): Reduced 
flow velocity along the average flapper position, Ub\ (b): Reduced kinetic energy, E. For both cases: r] a /rj = 0.1 (left solid line, 
blue); r/s/r] = 0.01 (red dot-dashed line); n a /rj = 0.001 (orange dotted line); r) a /ri = 0.0001 (right solid line, green); r\ s jr\ = 
(black dashed line). 



2. Kinetic energy 

Another possible optimization measure is related to the total kinetic energy of the average flow. Since the velocity 
field is directly proportional to the radius, it takes the general form (ui r ) = rf(9) and {uab) — r g(9), where the 
functions f(9) and g{9) can be found from Eq. (41). Therefore, the dependence of the total kinetic energy of the 
average flow upon the Deborah number can be characterized by a reduced energy given by the integral over the polar 
angle 

J3(De,C)= r /2 [f(9)] 2 + [g(e)] 2 d9, (44) 
Jo 

and is given analytically by 

De 6 7r(-i + C) 2 (i + C + 2De 2 C) 2 



£(De,C) 



98304 (1 + De 2 ) 4 (l + 4Dc 2 ) 2 (l + De 2 ( 2 ) 2 
-2074 - 4028Dc 2 - 2058Dc 4 + 1054tt 2 



+ 858DcV - 144De 4 7r 2 - 174tt 4 - 45De 2 7r 4 
+ 36Dc 4 tt 4 + 9tt 6 - 120Dc 2 C + 88De 4 C 
+ 1250DeVC + 1146De 4 7r 2 C - 303De 2 7r 4 C 
- 117De 4 7r 4 C + 18De 2 7r 6 C - 104De 4 C 2 

+ 52Dc 4 tt 2 C 2 - 93De 4 7r 4 C 2 + 9DeVc 2 ] . (45) 

With a fixed relative viscosity, for small values of De, we have E ~ De 6 , whereas E ~ 1/De 2 for large values of Dc, 
so an optimal De should exist. The function E(De, £) is plotted for different values of the relative viscosity in Fig. 5b, 
and similarly to the previous section we see indeed the existence of an optimal value of De for each £. 



3. Enstrophy 

Finally, we also consider the dependence of the enstrophy of the flow upon the Deborah number. The total enstrophy 
of the flow is proportional to the integral 

«r/2 

£(De,C) = / ojld9, (46) 
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FIG. 6: Optimal Deborah number, De opt , as a function of relative viscosity, rj s /r], by optimizing (a) the flow along the average 
position of the flapper 8 = n/2 (blue solid line) and (b) the flow kinetic energy (red dots). 



which can be analytically calculated to be 

Dc 6 7r(-l + C) 2 ( 1 + C + 2Dc 2 C) ; 



98304 (1 + De 2 ) 4 (l + 4De 2 ) 2 (l + De 2 C 2 ) 2 

1800 - 14256Dc 2 - 12744De 4 + 616tt 2 
+ 2424DcV + 1440DcV - 120tt 4 - 72Dc 2 tt 4 
+ 9tt 6 + 10656De 2 C + 11232De 4 C - 1192DeVC 
- 456De 4 7r 2 C - 168DeVC - 72De 4 7r 4 C 



18Dc 2 tt 6 C - 288De 4 C z - 368Dc 4 tt z C 



4X2 



- 48Dc 4 7r 4 C 2 + 9De 4 7r 6 C 2 



(47) 



The variation of £ with De turns out to be very similar to the one for E (Eq. 45, shown in Fig. 5b), and is not 
reproduced here. 



4- Optimal Deborah number 

We next compute numerically the optimal Deborah number, maximizing the pumping measures in both Eqs. (43) 
and (45), as a function of the relative viscosity. The results are displayed in Fig. 6. The optimality conditions turn out 
to quantitatively agree for both pumping measures, and correspond to an inverse linear relationship between De op t 
and C- This scaling is confirmed by an asymptotic analysis of the exact analytical formula for De op t found by setting 
the partial derivative of Eq. (43) to zero, and showing that indeed De opt ~ 1/C for small values of (. 



VI. DISCUSSION 



In this paper, we have considered what is arguably the simplest geometrical setup to demonstrate that net fluid 
pumping can be obtained from the purely sinusoidal forcing of a viscoelastic fluid. The fluid was modeled as an 
Oldroyd-B fluid both for simplicity and because of the physical relevance of the model. The main result we obtained 
is the time-averaged flow, described by Eq. (41), generated by the reciprocally flapping motion. In accordance with 
the scallop theorem, setting the Deborah number to zero in Eq. (41) leads to no flow, but a net flow occurs for all 
nonzero values of De. Our calculations allow us to demonstrate explicitly the breaking of the scallop theorem in the 
context of fluid pumping, and suggest the possibility of taking advantage of the intrinsic nonlinearities of complex 
fluids for their transport. Physically, such flow is being driven by normal-stress differences arising in the fluid and due 
to the stretching of the polymeric microstructures by the background flow. The calculation was done asymptotically 
for small-amplitude flapping, and the net flow occurs at fourth order. As in the classic work by Moffatt [49], our 
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results should be understood as similarity solutions which are valid close enough to the fixed hinge point such that 
the inertial effects are negligible. The advantage of such theoretical treatment is that it allows us to obtain the entire 
flow field analytically, in particular the spatial structure of the flow, and the dependancc of the net pumping on the 
actuation parameters (the flapping frequency) and the material properties of the fluid (relaxation time and viscosities). 
Taking advantage of these analytical results, we have been able to analytically optimize the pumping performance, 
and derive the optimal Deborah number as a function of the fluid ratio of solvent to total viscosity. Although we have 
considered here the simplest geometrical and dynamical setup possible, the results motivate future work which will 
focus on the flapping of three-dimensional finite-size appendages in polymeric fluids. 

We now turn to the relevance of our results to biological transport. In Newtonian fluids, only the non-reciprocal 
component of the motion of cilia - i.e. the difference between their effective and recovery strokes - affects fluid 
transport [5, 50]. In contrast, we show in this paper that the back- and- forth components of cilia motion, which is 
reciprocal, does influence transport in the case of viscoelastic biological fluids. The effect is expected to be crucial 
since the typical Deborah number in ciliary transport is large, and elastic effects of the fluid are therefore likely to be 
significant. For example, from rheological measurements [51-53], we know the relaxation time of respiratory mucus 
ranges between A w 30 — 100 s, and that of the cervical mucus present in female reproductive tract ranges from 
Awl — 100 s [51, 54]. In addition, cilia typically oscillate at frequencies of / = uj/2tt « 5 — 50 Hz [5], and therefore, 
ciliary transport of mucus occurs at large (or very large) Deborah numbers, De = \u> ~ 10 to 10 . 

In addition, the results of our paper should be contrasted with previous work. It was shown in Rcf. [21] that the 
presence of polymeric stresses leads to a decrease of the speed at which a fluid is pumped by a waving sheet - in that 
case a complex fluid led therefore to a degradation of the transport performance. In contrast, wc demonstrate in the 
current paper a mode of actuation which is rendered effective by the presence of polymeric stresses - the complex 
fluid leads therefore in this case to an improvement of the transport performance. For a general actuation gait, it is 
therefore not known a priori whether the presence of a complex fluid will lead to a degradation or an improvement 
of the pumping performance, and whether or not a general classification depending on the type of actuation gait can 
be derived remains a question to be addressed in the future. 
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